5 Practical Issues in the Specification and Estimation of Discrete Choice Models

library(knitr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(discrtr) # A companion package for the book Introduction to Discrete Choice Analysis with `R`
library(dplyr) # A Grammar of Data Manipulation
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
library(mlogit) # Multinomial Logit Models
## Cargando paquete requerido: dfidx
## 
## Adjuntando el paquete: 'dfidx'
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(readr) # Read Rectangular Text Data
#library(stargazer) # Well-Formatted Regression and Summary Statistics Tables
library(gplots) # Various R Programming Tools for Plotting Data
## 
## Adjuntando el paquete: 'gplots'
## 
## The following object is masked from 'package:stats':
## 
##     lowess
library(ggmosaic) # Mosaic Plots in the 'ggplot2' Framework
library(treemapify)
library(ggridges)
library(ggalluvial)
library(evd)
library(htmlwidgets) # HTML Widgets for R
library(kableExtra) # Construct Complex Table with kable and Pipe Syntax
## 
## Adjuntando el paquete: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(plotly) # Create interactive web graphics
## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(tidyr) # Tidy messy data
#library(webshot2) # Take screenshots of web pages

Funciones de utilidad

data("mc_commute_wide", package = "discrtr")
mc_commute_wide[1:6, 1:10]
##          id choice available.Cycle available.Walk available.HSR available.Car
## 1 566872636    HSR              No            Yes           Yes            No
## 2 566873140    HSR              No            Yes           Yes           Yes
## 3 566874266    HSR              No             No           Yes           Yes
## 4 566874842   Walk             Yes            Yes           Yes            No
## 5 566881170   Walk             Yes            Yes           Yes            No
## 6 566907438   Walk              No            Yes           Yes            No
##   time.Cycle time.Walk time.HSR time.Car
## 1         NA  21.31439        5       NA
## 2         NA  12.78863       10        2
## 3         NA        NA       15        4
## 4   5.828157  20.00000        8       NA
## 5   5.828157  20.00000        5       NA
## 6         NA  10.00000        3       NA

Se toma a continuacion choice como columna y tiempo las variables que se pueden elegir

example_wide <- mc_commute_wide %>% dplyr::select(id, choice, starts_with("time")) |> 
                filter(id %in% c(566910139, 566873140, 566872636))
example_wide
##          id choice time.Cycle time.Walk time.HSR time.Car
## 1 566872636    HSR         NA  21.31439        5       NA
## 2 566873140    HSR         NA  12.78863       10        2
## 3 566910139   Walk   4.371118  15.00000       20        5

Usaremos la misma selección de variables y filas para ver la diferencia con una tabla en formato por longitud.

example_wide |>  pivot_longer(cols = starts_with("time."), names_prefix = "time.", names_to="alternative", values_to="time")
## # A tibble: 12 × 4
##           id choice alternative  time
##        <dbl> <fct>  <chr>       <dbl>
##  1 566872636 HSR    Cycle       NA   
##  2 566872636 HSR    Walk        21.3 
##  3 566872636 HSR    HSR          5   
##  4 566872636 HSR    Car         NA   
##  5 566873140 HSR    Cycle       NA   
##  6 566873140 HSR    Walk        12.8 
##  7 566873140 HSR    HSR         10   
##  8 566873140 HSR    Car          2   
##  9 566910139 Walk   Cycle        4.37
## 10 566910139 Walk   Walk        15   
## 11 566910139 Walk   HSR         20   
## 12 566910139 Walk   Car          5

Nótese como el mismo id se repite en cuatro filas, pero el tiempo está dado en una sola. El package mlogit utiliza tablas en formato por longitud. Dado que las tablas en formato por amplitud son más comunes, este package incluye una función para cambiar la tabla. Para ello es necesario indicar cuáles son las variables que varían en cada alternativa, en este caso son cuatro variables: time, access (el tiempo necesario para llegar a una parada de autobús), wait (tiempo de espera por un autobús) y transfer (el número de transferencias cuando se utiliza el autobús). Las últimas tres variables son específicas de HSR (cuando se utiliza el transporte público)

example_long <- mc_commute_wide %>% filter(id %in% c(566910139, 566873140, 566872636)) |> 
                mlogit.data(shape="wide", choice="choice", varying=3:22)

example_long
## ~~~~~~~
##  first 10 observations out of 12 
## ~~~~~~~
##           id choice parking vehind gender age                          shared
## 1  566872636  FALSE      No     No    Man  21 Living in Shared Accommodations
## 2  566872636  FALSE      No     No    Man  21 Living in Shared Accommodations
## 3  566872636   TRUE      No     No    Man  21 Living in Shared Accommodations
## 4  566872636  FALSE      No     No    Man  21 Living in Shared Accommodations
## 5  566873140  FALSE      No    Yes    Man  23                              No
## 6  566873140  FALSE      No    Yes    Man  23                              No
## 7  566873140   TRUE      No    Yes    Man  23                              No
## 8  566873140  FALSE      No    Yes    Man  23                              No
## 9  566910139  FALSE      No    Yes    Man  49                              No
## 10 566910139  FALSE      No    Yes    Man  49                              No
##                family child street_density sidewalk_density      LAT      LONG
## 1                  No    No       14.37621         22.63322 43.26302 -79.90074
## 2                  No    No       14.37621         22.63322 43.26302 -79.90074
## 3                  No    No       14.37621         22.63322 43.26302 -79.90074
## 4                  No    No       14.37621         22.63322 43.26302 -79.90074
## 5                  No    No       19.49754         39.64003 43.25885 -79.90476
## 6                  No    No       19.49754         39.64003 43.25885 -79.90476
## 7                  No    No       19.49754         39.64003 43.25885 -79.90476
## 8                  No    No       19.49754         39.64003 43.25885 -79.90476
## 9  Living with Family   Yes       15.39380         26.06793 43.25597 -79.91260
## 10 Living with Family   Yes       15.39380         26.06793 43.25597 -79.91260
##    PersonalVehComf_SD PersonalVehComf_D PersonalVehComf_A PersonalVehComf_SA
## 1                   0                 0                 1                  0
## 2                   0                 0                 1                  0
## 3                   0                 0                 1                  0
## 4                   0                 0                 1                  0
## 5                   0                 0                 0                  0
## 6                   0                 0                 0                  0
## 7                   0                 0                 0                  0
## 8                   0                 0                 0                  0
## 9                   0                 0                 1                  0
## 10                  0                 0                 1                  0
##    Fun_SD Fun_D Fun_A Fun_SA ActiveNeigh_SD ActiveNeigh_D ActiveNeigh_A
## 1       0     1     0      0              0             0             1
## 2       0     1     0      0              0             0             1
## 3       0     1     0      0              0             0             1
## 4       0     1     0      0              0             0             1
## 5       0     0     0      1              0             0             0
## 6       0     0     0      1              0             0             0
## 7       0     0     0      1              0             0             0
## 8       0     0     0      1              0             0             0
## 9       0     0     1      0              0             0             0
## 10      0     0     1      0              0             0             0
##    ActiveNeigh_SA UsefulTrans_SD UsefulTrans_D UsefulTrans_A UsefulTrans_SA
## 1               0              0             0             1              0
## 2               0              0             0             1              0
## 3               0              0             0             1              0
## 4               0              0             0             1              0
## 5               0              1             0             0              0
## 6               0              1             0             0              0
## 7               0              1             0             0              0
## 8               0              1             0             0              0
## 9               0              0             0             1              0
## 10              0              0             0             1              0
##    BusComf_SD BusComf_D BusComf_A BusComf_SA TravelAlone_SD TravelAlone_D
## 1           0         1         0          0              1             0
## 2           0         1         0          0              1             0
## 3           0         1         0          0              1             0
## 4           0         1         0          0              1             0
## 5           0         0         0          0              1             0
## 6           0         0         0          0              1             0
## 7           0         0         0          0              1             0
## 8           0         0         0          0              1             0
## 9           0         0         1          0              0             1
## 10          0         0         1          0              0             1
##    TravelAlone_A TravelAlone_SA Shelters_SD Shelters_D Shelters_A Shelters_SA
## 1              0              0           0          0          1           0
## 2              0              0           0          0          1           0
## 3              0              0           0          0          1           0
## 4              0              0           0          0          1           0
## 5              0              0           0          0          1           0
## 6              0              0           0          0          1           0
## 7              0              0           0          0          1           0
## 8              0              0           0          0          1           0
## 9              0              0           0          0          0           0
## 10             0              0           0          0          0           0
##    Community_SD Community_D Community_A Community_SA personal_veh_comfortable
## 1             0           0           0            0                        A
## 2             0           0           0            0                        A
## 3             0           0           0            0                        A
## 4             0           0           0            0                        A
## 5             0           1           0            0                        N
## 6             0           1           0            0                        N
## 7             0           1           0            0                        N
## 8             0           1           0            0                        N
## 9             0           0           1            0                        A
## 10            0           0           1            0                        A
##    getting_there_fun like_active_neighborhood commute_useful_transition
## 1                  D                        A                         A
## 2                  D                        A                         A
## 3                  D                        A                         A
## 4                  D                        A                         A
## 5                 SA                        N                        SD
## 6                 SA                        N                        SD
## 7                 SA                        N                        SD
## 8                 SA                        N                        SD
## 9                  A                        N                         A
## 10                 A                        N                         A
##    buses_comfortable prefer_travel_alone shelter_good_quality sense_community
## 1                  D                  SD                    A               N
## 2                  D                  SD                    A               N
## 3                  D                  SD                    A               N
## 4                  D                  SD                    A               N
## 5                  N                  SD                    A               D
## 6                  N                  SD                    A               D
## 7                  N                  SD                    A               D
## 8                  N                  SD                    A               D
## 9                  A                   D                    N               A
## 10                 A                   D                    N               A
##    numna   alt available      time access wait transfer chid    idx
## 1      2   Car        No        NA      0    0        0    1  1:Car
## 2      2 Cycle        No        NA      0    0        0    1 1:ycle
## 3      2   HSR       Yes  5.000000      3   15        0    1  1:HSR
## 4      2  Walk       Yes 21.314387      0    0        0    1 1:Walk
## 5      3   Car       Yes  2.000000      0    0        0    2  2:Car
## 6      3 Cycle        No        NA      0    0        0    2 2:ycle
## 7      3   HSR       Yes 10.000000      4   15        0    2  2:HSR
## 8      3  Walk       Yes 12.788632      0    0        0    2 2:Walk
## 9      4   Car       Yes  5.000000      0    0        0    3  3:Car
## 10     4 Cycle       Yes  4.371118      0    0        0    3 3:ycle
## 
## ~~~ indexes ~~~~
##    chid   alt
## 1     1   Car
## 2     1 Cycle
## 3     1   HSR
## 4     1  Walk
## 5     2   Car
## 6     2 Cycle
## 7     2   HSR
## 8     2  Walk
## 9     3   Car
## 10    3 Cycle
## indexes:  1, 2
data.frame(example_long) |>  dplyr::select(id, choice, alt, starts_with("time"), idx)
##           id choice   alt      time    idx
## 1  566872636  FALSE   Car        NA  1:Car
## 2  566872636  FALSE Cycle        NA 1:ycle
## 3  566872636   TRUE   HSR  5.000000  1:HSR
## 4  566872636  FALSE  Walk 21.314387 1:Walk
## 5  566873140  FALSE   Car  2.000000  2:Car
## 6  566873140  FALSE Cycle        NA 2:ycle
## 7  566873140   TRUE   HSR 10.000000  2:HSR
## 8  566873140  FALSE  Walk 12.788632 2:Walk
## 9  566910139  FALSE   Car  5.000000  3:Car
## 10 566910139  FALSE Cycle  4.371118 3:ycle
## 11 566910139  FALSE   HSR 20.000000  3:HSR
## 12 566910139   TRUE  Walk 15.000000 3:Walk

Dado que hay cuatro alternativas, cada fila corresponde a la situación de elección para una alternativa para un individuo (las tablas en formato por longitud por lo general tienen filas por cada decisor). El conjunto de índices contiene dos variables, la primera identifica al individuo y la segunda la alternativa:

data.frame(example_long$idx)
##    chid   alt
## 1     1   Car
## 2     1 Cycle
## 3     1   HSR
## 4     1  Walk
## 5     2   Car
## 6     2 Cycle
## 7     2   HSR
## 8     2  Walk
## 9     3   Car
## 10    3 Cycle
## 11    3   HSR
## 12    3  Walk

En este caso las variables disponibles para el análisis son:

colnames(example_long)
##  [1] "id"                        "choice"                   
##  [3] "parking"                   "vehind"                   
##  [5] "gender"                    "age"                      
##  [7] "shared"                    "family"                   
##  [9] "child"                     "street_density"           
## [11] "sidewalk_density"          "LAT"                      
## [13] "LONG"                      "PersonalVehComf_SD"       
## [15] "PersonalVehComf_D"         "PersonalVehComf_A"        
## [17] "PersonalVehComf_SA"        "Fun_SD"                   
## [19] "Fun_D"                     "Fun_A"                    
## [21] "Fun_SA"                    "ActiveNeigh_SD"           
## [23] "ActiveNeigh_D"             "ActiveNeigh_A"            
## [25] "ActiveNeigh_SA"            "UsefulTrans_SD"           
## [27] "UsefulTrans_D"             "UsefulTrans_A"            
## [29] "UsefulTrans_SA"            "BusComf_SD"               
## [31] "BusComf_D"                 "BusComf_A"                
## [33] "BusComf_SA"                "TravelAlone_SD"           
## [35] "TravelAlone_D"             "TravelAlone_A"            
## [37] "TravelAlone_SA"            "Shelters_SD"              
## [39] "Shelters_D"                "Shelters_A"               
## [41] "Shelters_SA"               "Community_SD"             
## [43] "Community_D"               "Community_A"              
## [45] "Community_SA"              "personal_veh_comfortable" 
## [47] "getting_there_fun"         "like_active_neighborhood" 
## [49] "commute_useful_transition" "buses_comfortable"        
## [51] "prefer_travel_alone"       "shelter_good_quality"     
## [53] "sense_community"           "numna"                    
## [55] "alt"                       "available"                
## [57] "time"                      "access"                   
## [59] "wait"                      "transfer"                 
## [61] "chid"                      "idx"

Iniciaremos definiendo una fórmula que sólo considera el tiempo de transporte.

# Function `mFormula()` is used to define multi-part formulas of the form:
# y ~ x | z | w, which in the notation used for the anatomy of utility functions is
# choice ~ alternative vars. with generic coefficients |
# individual vars. with specific coefficients |
# alternative vars. with specific coefficients
# In this formula time is one of x variables
f1 <- mFormula(choice ~ time)

La función model.matrix puede ser usada para examinar cómo la fórmula es aplicada a los datos:

f1 |> model.matrix(example_long)
##   (Intercept):Cycle (Intercept):HSR (Intercept):Walk      time
## 1                 0               1                0  5.000000
## 2                 0               0                1 21.314387
## 3                 0               0                0  2.000000
## 4                 0               1                0 10.000000
## 5                 0               0                1 12.788632
## 6                 0               0                0  5.000000
## 7                 1               0                0  4.371118
## 8                 0               1                0 20.000000
## 9                 0               0                1 15.000000

Definamos una fórmula con una variable específica de cada individuo

# Function `mFormula()` is used to define multi-part formulas of the form:
# y ~ x | z | w, which in the notation used for the anatomy of utility functions is
# choice ~ alternative vars. with generic coefficients |
# individual vars. with specific coefficients |
# alternative vars. with specific coefficients
# In this formula `time` is one of x variables and `sidewalk_density` is one of z variables
f2 <- mFormula(choice ~ time | sidewalk_density)

La matriz del modelo está dada por el tiempo que se toma usar lod diferentes medios de transporte:

f2 |>  model.matrix(example_long)
##   (Intercept):Cycle (Intercept):HSR (Intercept):Walk      time
## 1                 0               1                0  5.000000
## 2                 0               0                1 21.314387
## 3                 0               0                0  2.000000
## 4                 0               1                0 10.000000
## 5                 0               0                1 12.788632
## 6                 0               0                0  5.000000
## 7                 1               0                0  4.371118
## 8                 0               1                0 20.000000
## 9                 0               0                1 15.000000
##   sidewalk_density:Cycle sidewalk_density:HSR sidewalk_density:Walk
## 1                0.00000             22.63322               0.00000
## 2                0.00000              0.00000              22.63322
## 3                0.00000              0.00000               0.00000
## 4                0.00000             39.64003               0.00000
## 5                0.00000              0.00000              39.64003
## 6                0.00000              0.00000               0.00000
## 7               26.06793              0.00000               0.00000
## 8                0.00000             26.06793               0.00000
## 9                0.00000              0.00000              26.06793

Ahora consideramos que el tiempo es específico de cada alternativa, en lugar de tener coeficientes genéricos, por medio de la fórmula:

f3 <- mFormula(choice ~ 0 | sidewalk_density | time)

Nótese que, dado que no definimos variables específicas de cada alternativa con coeficientes genéricos, tenemos que establecerlas como cero. La matriz del modelo es:

f3 %>% model.matrix(example_long)
##   (Intercept):Cycle (Intercept):HSR (Intercept):Walk sidewalk_density:Cycle
## 1                 0               1                0                0.00000
## 2                 0               0                1                0.00000
## 3                 0               0                0                0.00000
## 4                 0               1                0                0.00000
## 5                 0               0                1                0.00000
## 6                 0               0                0                0.00000
## 7                 1               0                0               26.06793
## 8                 0               1                0                0.00000
## 9                 0               0                1                0.00000
##   sidewalk_density:HSR sidewalk_density:Walk time:Car time:Cycle time:HSR
## 1             22.63322               0.00000        0   0.000000        5
## 2              0.00000              22.63322        0   0.000000        0
## 3              0.00000               0.00000        2   0.000000        0
## 4             39.64003               0.00000        0   0.000000       10
## 5              0.00000              39.64003        0   0.000000        0
## 6              0.00000               0.00000        5   0.000000        0
## 7              0.00000               0.00000        0   4.371118        0
## 8             26.06793               0.00000        0   0.000000       20
## 9              0.00000              26.06793        0   0.000000        0
##   time:Walk
## 1   0.00000
## 2  21.31439
## 3   0.00000
## 4   0.00000
## 5  12.78863
## 6   0.00000
## 7   0.00000
## 8   0.00000
## 9  15.00000

Experimento 1 Consiste en que los valores mu y beta son iguales

ts <- tibble(Individual=1:6, Choice=c("A","A", "B", "A", "B", "B"),
             yA=c(1,1,0,1,0,0), yB=c(0, 0, 1, 0, 1, 1),
             xA=c(5,2,5,1,4,3), xB=c(4, 5, 2, 6, 1, 4))


# Set the parameters:
mu <- 0
beta <- 0

my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))  
}

# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A

# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A

# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A

# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A

# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A

# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A

# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]

# Create data frame to tabulate results:
df_experiment_1 <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))
# Display table
kable(df_experiment_1, "html", digits = 4, booktabs = TRUE, align = c("c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
footnote(general = paste("The value of the likelihood function in Example 1 is: ",
round(L, digits = 4)))
Individual Choice PA PB
1 A 0.5 0.5
2 A 0.5 0.5
3 B 0.5 0.5
4 A 0.5 0.5
5 B 0.5 0.5
6 B 0.5 0.5
Note:
The value of the likelihood function in Example 1 is: 0.0156

Experimento 2 1 Consiste en que los valores mu y beta son diferentes

# Set the parameters:
mu <- 0.5
beta <- -0.5

my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))  
}

# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A

# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A

# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A

# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A

# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A

# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A

# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]

# Create data frame to tabulate results:
df_experiment_2 <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))
# Display table
kable(df_experiment_2, "html", digits = 4, booktabs = TRUE, align = c("c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
footnote(general = paste("The value of the likelihood function in Example 2 is: ",
round(L, digits = 4)))
Individual Choice PA PB
1 A 0.2689 0.7311
2 A 0.7311 0.2689
3 B 0.1192 0.8808
4 A 0.8808 0.1192
5 B 0.1192 0.8808
6 B 0.5000 0.5000
Note:
The value of the likelihood function in Example 2 is: 0.0672

Al cambiar los valores provoca que cambien las probabilidades

# Create a grid to plot the likelihood function
mu = seq(from = -1, to = 1, by = 0.01)
beta = seq(from = -2, to = 0, by = 0.01)
coeffs <- expand.grid(mu, beta)

my_prob <- function(xA, xB, mu, beta){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))  
}


ts <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
       Choice = c("A", "A", "B", "A", "B", "B"),
       yA = c(1, 1, 0, 1, 0, 0),
       yB = c(0, 0, 1, 0, 1, 1),
       xA = c(5, 2, 5, 1, 4, 3),
       xB = c(4, 5, 2, 6, 1, 4))


# Define the likelihood function
lkh <- function(mu = 0, beta = 0){
       
       P1A <- my_prob(ts$xA[1], ts$xB[1], mu, beta)
       P1B <- 1 - P1A

       P2A <- my_prob(ts$xA[2], ts$xB[2], mu, beta)
       P2B <- 1 - P2A

       P3A <- my_prob(ts$xA[3], ts$xB[3], mu, beta)
       P3B <- 1 - P3A

       P4A <- my_prob(ts$xA[4], ts$xB[4], mu, beta)
       P4B <- 1 - P4A

       P5A <- my_prob(ts$xA[5], ts$xB[5], mu, beta)
       P5B <- 1 - P5A

       P6A <- my_prob(ts$xA[6], ts$xB[6], mu, beta)
       P6B <- 1 - P6A

       P1A^ts$yA[1] * P1B^ts$yB[1] *
       P2A^ts$yA[2] * P2B^ts$yB[2] *
       P3A^ts$yA[3] * P3B^ts$yB[3] *
       P4A^ts$yA[4] * P4B^ts$yB[4] *
       P5A^ts$yA[5] * P5B^ts$yB[5] *
       P6A^ts$yA[6] * P6B^ts$yB[6]
}

# Evaluate the likelihood function on the grid
L <- lkh(mu = coeffs$Var1, beta = coeffs$Var2)
L <- data.frame(mu = coeffs$Var1, beta = coeffs$Var2, L)
L <- xtabs(L ~ beta + mu, L) %>% # Convert to cross-tabulation matrix
unclass() # Drop the xtabs class (plotly does not like it)

likelihood_plot <- plot_ly(z = ~L, x = ~mu, y = ~beta) %>%
add_surface() %>%
layout(scene = list(
xaxis = list(title = "x-axis (mu)"),
yaxis = list(title = "y-axis (beta)"),
zaxis = list(title = "z-axis (L)")))

likelihood_plot

En la gráfica se puede observar que los valores aproximados que maximizan la función de verosimilitud son mu = 0.10 y bet = 0.65. Si utilizamos estos coeficientes para calcular las probabilidades, podemos comparar las probabilidades de los experimentos 1 y 2:

# Set the parameters:
mu <- 0.1
beta <- -0.65

my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))  
}

# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A

# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A

# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A

# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A

# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A

# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A

# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]

# Create data frame to tabulate results:
df_approx_solution <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))

# Join tables for displaying results
df <- df_experiment_1 %>% left_join(df_experiment_2,
      by = c("Individual", "Choice")) %>%
      left_join(df_approx_solution, 
      by = c("Individual", "Choice"))



# Display table
kable(df,"html", digits = 4, booktabs = TRUE, col.names = c("Individual", "Choice", "PA", "PB", "PA", "PB", "PA", "PB"),
      align = c("l", "c", "c", "c", "c", "c", "c", "c")) %>%
      kable_styling(latex_options = c("striped")) %>%
      add_header_above(c(" " = 1, " " = 1,
      "Experiment 1" = 2,
      "Experiment 2" = 2,
      "Approx Max Likelihood" = 2))%>%
footnote(general = paste("The approximate optimal value of the likelihood function is: ",
round(L, digits = 4)))
Experiment 1
Experiment 2
Approx Max Likelihood
Individual Choice PA PB PA PB PA PB
1 A 0.5 0.5 0.2689 0.7311 0.3208 0.6792
2 A 0.5 0.5 0.7311 0.2689 0.8641 0.1359
3 B 0.5 0.5 0.1192 0.8808 0.1141 0.8859
4 A 0.5 0.5 0.8808 0.1192 0.9589 0.0411
5 B 0.5 0.5 0.1192 0.8808 0.1141 0.8859
6 B 0.5 0.5 0.5000 0.5000 0.6341 0.3659
Note:
The approximate optimal value of the likelihood function is: 0.0763

Dado que la función de verosimilitud toma valores entre 0 y 1,la función log-likelihood toma valores menores que cero.

modelo logit

Consideraremos nuevamente el modelo para elección de modo de transporte. Es necesario transformar todo el conjunto de datos del formato por amplitud por longitud:

mc_commute_long <- mc_commute_wide |> 
mlogit.data(shape="wide",
            # Name of column with the choices
            choice = "choice",
            # Numbers of columns with attributes that vary by alternative
            varying = 3:22)

Esta función requiere al menos dos argumentos: una mFormula y un conjunto de datos. Podemos verificar que las fórmulas definidas previamente son objetos de este tipo:

class(f1)
## [1] "mFormula" "Formula"  "formula"
class(f2)
## [1] "mFormula" "Formula"  "formula"
class(f3)
## [1] "mFormula" "Formula"  "formula"
# Function `mlogit()` is used to estimate logit models
# It needs a multi-part formula and a data set in long form
model1 <- mlogit(f1, mc_commute_long)
# Function `summary()` give the summary of data objects,
# including the output of model estimation algorithms
summary(model1)
## 
## Call:
## mlogit(formula = choice ~ time, data = mc_commute_long, method = "nr")
## 
## Frequencies of alternatives:choice
##      Car    Cycle      HSR     Walk 
## 0.204364 0.034182 0.244364 0.517091 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000102 
## successive function values within tolerance limits 
## 
## Coefficients :
##                    Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept):Cycle  0.366898   0.201434  1.8214   0.06854 .  
## (Intercept):HSR    0.707267   0.118104  5.9885 2.118e-09 ***
## (Intercept):Walk   3.212834   0.184169 17.4451 < 2.2e-16 ***
## time              -0.056494   0.005685 -9.9374 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -768.63
## McFadden R^2:  0.50323 
## Likelihood ratio test : chisq = 1557.2 (p.value = < 2.22e-16)

Los resultados incluyen la proporción de frecuencia observada de cada alternativa además de información acerca del procedimiento de optimización.

Ahora se estima otro modelo usando la segunda fórmula:

model2 <- mlogit(f2, mc_commute_long)

summary(model2)
## 
## Call:
## mlogit(formula = choice ~ time | sidewalk_density, data = mc_commute_long, 
##     method = "nr")
## 
## Frequencies of alternatives:choice
##      Car    Cycle      HSR     Walk 
## 0.204364 0.034182 0.244364 0.517091 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000287 
## successive function values within tolerance limits 
## 
## Coefficients :
##                          Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept):Cycle      -0.3058620  0.4541379 -0.6735  0.500629    
## (Intercept):HSR         0.6761001  0.2150745  3.1436  0.001669 ** 
## (Intercept):Walk        2.3453233  0.3124706  7.5057 6.106e-14 ***
## time                   -0.0562191  0.0056910 -9.8785 < 2.2e-16 ***
## sidewalk_density:Cycle  0.0273085  0.0169331  1.6127  0.106803    
## sidewalk_density:HSR    0.0016950  0.0085746  0.1977  0.843299    
## sidewalk_density:Walk   0.0350392  0.0111749  3.1355  0.001715 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -760.42
## McFadden R^2:  0.50854 
## Likelihood ratio test : chisq = 1573.7 (p.value = < 2.22e-16)

Es posible seleccionar el nivel de referencia para las utilidades cuando se estima el modelo, en el siguiente ejemplo se utiliza la alternativa “walk” como referencia:

model2 <- mlogit(f2, mc_commute_long, reflevel = "Walk")

summary(model2)
## 
## Call:
## mlogit(formula = choice ~ time | sidewalk_density, data = mc_commute_long, 
##     reflevel = "Walk", method = "nr")
## 
## Frequencies of alternatives:choice
##     Walk      Car    Cycle      HSR 
## 0.517091 0.204364 0.034182 0.244364 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000287 
## successive function values within tolerance limits 
## 
## Coefficients :
##                          Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept):Car        -2.3453233  0.3124706 -7.5057 6.106e-14 ***
## (Intercept):Cycle      -2.6511853  0.4187557 -6.3311 2.434e-10 ***
## (Intercept):HSR        -1.6692232  0.2302925 -7.2483 4.221e-13 ***
## time                   -0.0562191  0.0056910 -9.8785 < 2.2e-16 ***
## sidewalk_density:Car   -0.0350392  0.0111749 -3.1355  0.001715 ** 
## sidewalk_density:Cycle -0.0077307  0.0148205 -0.5216  0.601935    
## sidewalk_density:HSR   -0.0333442  0.0084297 -3.9556 7.635e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -760.42
## McFadden R^2:  0.50854 
## Likelihood ratio test : chisq = 1573.7 (p.value = < 2.22e-16)

Observemos que ahora dos coeficientes de sidewalk son significativos. Mientras la densidad de acera no cambia significativamente la utilidad entre utilizar bicicleta o caminar, vivir en un lugar con alta densidad de aceras reduce la utilidad de usar automóvil y transporte público respecto a caminar.

summary(mc_commute_long$sidewalk_density)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   18.19   22.63   24.18   35.70   59.41

Posteriormente, del conjunto de datos utilizados para estimar el modelo seleccionamos algunas observaciones para explorar las densidades de aceras en el rango entre 0 y 60, en intervalos de longitud 5. Por lo tanto, se seleccionan 52 observaciones (filas) de la tabla (trece niveles de densidad de aceras por cuatro alternativas):

mc_commute_predict <- mc_commute_long[1:52, ]

Reemplazamos la variable sidewalk density usando valores entre 0 y 60, en intervalos de longitud 5. Dado que cada alternativa es una fila, es necesario crear una sucesión de valores repetidos de la siguiente manera:

mc_commute_predict$sidewalk_density <- rep(seq(from=0, to=60, by=5), each=4)

Podemos examinar los valores asignados para sidewalk density

mc_commute_predict |>  data.frame() |>  select(sidewalk_density) |>  slice_head(n=10)
##    sidewalk_density
## 1                 0
## 2                 0
## 3                 0
## 4                 0
## 5                 5
## 6                 5
## 7                 5
## 8                 5
## 9                10
## 10               10

El conjunto de datos para predicción ahora incluye el rango de valores de sidewalk density en los cuales se tiene interés. El objetivo de la simulación es calcular la probabilidad de escoger alguna de las opciones de transporte conforme la densidad de acera varía, para ello también es necesario fijar los valores de otras variables, en particular se fija el valor del tiempo (se considera la mediana):

median(mc_commute_predict$time, na.rm=TRUE)
## [1] 10
mc_commute_predict$time <- 10

Ahora se calculan las probabilidades por medio de la función predict y el modelo 2:

probs <- predict(model2, newdata=mc_commute_predict)

El resultado (value) de predict es una matriz que contiene las probabilidades para trece niveles de densidad de acera y cuatro alternativas de transporte:

print(probs)
##         Walk        Car      Cycle        HSR
## 1  0.7381288 0.07072471 0.05208792 0.13905856
## 2  0.7646590 0.06149222 0.05191414 0.12193469
## 3  0.7887761 0.05323772 0.05152105 0.10646509
## 4  0.8105461 0.04591523 0.05093563 0.09260301
## 5  0.8300785 0.03946495 0.05018525 0.08027133
## 6  0.8475135 0.03381834 0.04929656 0.06937162
## 7  0.8630106 0.02890248 0.04829466 0.05979230
## 8  0.8767385 0.02464350 0.04720262 0.05141540
## 9  0.8888677 0.02096925 0.04604115 0.04412190
## 10 0.8995646 0.01781113 0.04482852 0.03779578
## 11 0.9089873 0.01510533 0.04358057 0.03232680
## 12 0.9172834 0.01279350 0.04231084 0.02761229
## 13 0.9245881 0.01082299 0.04103075 0.02355815

Para facilitar el procedimiento al realizar la gráfica, se agregan los valores de la densidad de acera y se cambia la forma de la tabla (por longitud), de tal manera que cada fila sea la probabilidad de la combinación algún modo de transporte y densidad de acera:

probs <- data.frame(sidewalk_density=seq(from=0, to=60, by=5), probs) %>% 
                    pivot_longer(cols=-sidewalk_density,
                                 names_to="Mode",
                                 values_to = "Probability")

probs %>% slice_head(n=10)
## # A tibble: 10 × 3
##    sidewalk_density Mode  Probability
##               <dbl> <chr>       <dbl>
##  1                0 Walk       0.738 
##  2                0 Car        0.0707
##  3                0 Cycle      0.0521
##  4                0 HSR        0.139 
##  5                5 Walk       0.765 
##  6                5 Car        0.0615
##  7                5 Cycle      0.0519
##  8                5 HSR        0.122 
##  9               10 Walk       0.789 
## 10               10 Car        0.0532

Ahora podemos generar la gráfica

ggplot(probs)+
  geom_line(aes(x=sidewalk_density, y=Probability, color=Mode), linewidth=1.2)+
  labs(x= expression("Sidewalk density (km/km"^2*")"),  y="Probability")

es denominado el modelo “Market Shares”:

f0 <- mFormula(choice ~ 1)
model0 <- mlogit(f0, mc_commute_long)
summary(model0)
## 
## Call:
## mlogit(formula = choice ~ 1, data = mc_commute_long, method = "nr")
## 
## Frequencies of alternatives:choice
##      Car    Cycle      HSR     Walk 
## 0.204364 0.034182 0.244364 0.517091 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 9.89E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##                    Estimate Std. Error  z-value Pr(>|z|)    
## (Intercept):Cycle -1.788207   0.157592 -11.3470  < 2e-16 ***
## (Intercept):HSR    0.178756   0.080839   2.2113  0.02702 *  
## (Intercept):Walk   0.928318   0.070464  13.1743  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1547.2
## McFadden R^2:  0 
## Likelihood ratio test : chisq = 0 (p.value = 1)

La verosimilitud logarítmica del modelo Market Shares es -1547.2, mientras que, para el modelo 2 es -760.418. Luego, la de McFadden para el modelo 2 es:

1-as.numeric(model2$logLik)/as.numeric(model0$logLik)
## [1] 0.5085353

En el resumen de los modelos, la prueba de proporción de verosimilitud está reportada, y se contrasta con el modelo nulo. Cuando es necesario comparar dos modelos no nulos se utiliza la función lrtest:

lrtest(model1, model2)
## Likelihood ratio test
## 
## Model 1: choice ~ time
## Model 2: choice ~ time | sidewalk_density
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   4 -768.63                         
## 2   7 -760.42  3 16.429  0.0009258 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1